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SUMMARY 


A realization of a stochastic atmosphere model for use in 
simulations is presented. The model provides pressure, density, 
temperature, and wind velocity as a function of latitude, longitude, and 
altitude, and is implemented in a three-degree of freedom simulation 
package. This implementation is used in the Monte Carlo simulation of an 
aeroassisted orbital transfer maneuver and results are compared to those 
of a more traditional approach. 


ABBREVIATIONS AND SYMBOLS 


Roman 

g 

h 

L 

P(n) 

R 

P 

T(n) 

T(n) 

u(n) 

v(n) 

z 


Gaussian random number 
altitude (km) 

correlation scale distance (km) 
pressure at time step n (N/m ) 
density autocorrelation 
temperature realized at time step n (K) 
mean temperature at time step n (K) 
eastward wind at step n (m/sec) 
northward wind at step n (m/sec) 
distance measure (km) 


Greek 

p(n) 

p L (n) 

p s (n) 

p(n) 

P* 

<r 

P 

<t> 


3 

realized density at time step n (kg/m ) 

large scale density perturbation at time step n (per 

mill) 

small scale density perturbation at time step n (per 
mill) 

mean density at time step n (kg/m ) 
total density perturbation p^ + (per mill) 
density standard deviation 
latitude (degrees) 


Subscripts 

H 


horizontal 


A /92-3/ys 7#- 



l large scale 

n time step n 

p pressure 

s small scale 

T temperature 

u wind speed 

p density 

k either l or s, as appropriate 

Miscellaneous 

(overbar) mean value 

’ (prime) total perturbation 


INTRODUCTION 

In order to perform rigorous simulations of the trajectories of 
advanced aerospace vehicles, it is desirable to have an accurate model 
of the Earth’s atmosphere. Unfortunately for the engineers designing 
such vehicles, the atmosphere is much too large and complex to 
completely determine its state at any given time. Such a system is best 
modeled as a stochastic process. 

This report documents an implementation of the Global Reference 
Atmosphere Model (GRAM) [1-2] as an option for the Program to Optimize 
Simulated Trajectories (POST) [3]. A brief description of the random 
atmosphere model is given in the first section. In the next section, 
implementation of the model in POST is discussed. The following section 
describes an application of the model to a Monte Carlo simulation. The 
final section presents conclusions. An appendix gives the detailed 
source code. 


ATMOSPHERIC SIMULATION 


The GRAM model was taken as the baseline for the atmospheric 
simulation development. The GRAM model describes the Earth’s atmosphere 
as a tabular function of altitude, longitude, latitude, and calendar 
month, plus a random component which accounts for such effects as 
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gravity waves, tides and turbulence. 

One aspect of the model which is important for trajectory 
simulations is the spatial correlation of the random variations. In the 
atmosphere, the pressure and density at a given point have values close 
to that of the pressure and density at a nearby point. Thus, in the 
model, some correlation function is needed to limit gradients in 
pressure and density to realistic values. 

Atmospheric properties at any point along a trajectory are given as 
a mean at that location plus random components which are correlated with 
the previous trajectory point [2]. For density, 

p(n) = P L (n) + p g (n) + p(n) (1) 

where p(n) is the density at the n th time step, p^ and p g are the large 
and small scale density perturbations, and p is the mean density. The 
large and small scale perturbations are so called because of their 
relative wavelengths. The large scale perturbations are tightly 
correlated so they vary gradually with spatial displacement, while the 
small scale perturbations are loosely correlated and vary more rapidly. 
In (1), 


p (n) = p(n)[R cr (n)(p(n-l)-p(n-l)/(<r (n-l)p(n-l)) 

K pK pR 

+ cr (n) g(l-R 2 ) 1/2 ] (2) 

pK pK 

where o' (n) is the standard deviation of density at the present point, 

pK 

o' (n-1) is its standard deviation from the previous time step, p(n-l) 

pK 

is the realized density from the previous time step, and g is a Gaussian 
random variable with zero mean and unity standard deviation. Throughout 
the document, the subscript k can be replaced consistently by l or s for 
large or small scale, as appropriate, 

Note that the o'* s for all parameters are stored in tables as a 
total variance, s, and a fraction of the variance due to large scale 
effects, f L » so 

°L = ^ f L S (3) 

o' = v (1 - f ) s (4) 

S L 
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An equation similar to (2) is used to determine P g (n). The large 
scale density autocorrelation used in (2) is 


R = exp[“ ((4x/L H )Z * (4Z/L » )2)1/2 1 

PK p K pK 


( 5 ) 


Here Ax is the horizontal distance in kilometers from the last 
trajectory point to the present point, Az is the vertical distance, and 
L and L are scale distances determined by altitude and latitude as 

H V 

PK pK 

discussed below. The correlations R and R are used to compute the 

tl ts r 

temperature variation, and are calculated in a similar manner. 

In the calculation of the random temperature 


T(n) = T(n) + T (n) + T (n) 
L s 

where 


T K (n) = T(n)[C K (T(n-l)-T(n-l))/T(n-l) 


with 


+ D (p(n)-p(n) )/p(n) + E g] 

K K 


c = 

K 


Ic (n)/o- (n-1 ) 3 [ (R - R R R )] 

T T TK pK 22 11 

(1-R 2 R 2 ) 

pK 11 


( 6 ) 


(7) 

( 8 ) 


d k = I ^^^V n - 1)1/( V»iV (n)) 


(9) 


E r = [o- 2 (n) - cV(n-l) - D 2 <r 2 (n) 

K T K T K p 


- 2C D R R it (n*l )<r (n) ] 

K K pK 11 T pK 


1/2 


( 10 ) 


where 

[(o' n (n)/P(n)) 2 - (o* (n)/p(n) ) 2 - (<r (n)/T(n)) 2 ] 

H = — S BL 15 (11) 

2(o- pK (n)/p(n) ) (o- TK (n)/T(n) ) 

and R has the same functional form as R , but is evaluated at step 
11 22 * 
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n-1. 

The pressure variation is determined from the temperature and 
density variations by expanding the ideal gas law to first order: 


P(n) = P(n) + P (n) + P (n) 
s L 


( 12 ) 


where 


P L (n) = P(n) [ (p L (n)/p(n) ) + (T L (n)/T(n) ) ] (13) 

and similarly for P g (n). 

The eastward wind component is given by: 


u = u + u + u 
l s 

with 


(14) 


u (n) = u(n) F u (n-1) + G p (n) + H g (15) 

K I K K K K K I 

where u(n) is the perturbed eastward wind at trajectory point n, u(n) is 
the mean east wind at point n, u^fn-l) is the east wind perturbation 
from the previous point, p R (n) is either the large or small scale 
density perturbation calculated from (2). The coefficients F,G, and H 
are given by: 


F = 


[<r,(n)/«r (n-l)][(R -R R R )] 
0 U UK pti U2p2 uipi 

(1-R 2 R 2 ) 
pK uipi 


(16) 


G = [R <r (n)-F <r (n-l)]/(R R <r (n)) 

K UK U K U PK Uipi pK 


(17) 


h k - i% 2(n » - f.V 1 "- 11 - G t V (n) 


- 2F G R R <r (n-l)cr (n)] 
K K PK uipi U PK 


1/2 


( 18 ) 


where R is given by a form analogous to (5) and R and R are 

UK U2p2 uipi 

taken from tables. 


The northward wind component is given by: 
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with 


v 


v + v + v 
L s 


(19) 


v (n) = v(n) I v (n-1 ) + J u (n) + Kg (20) 

I 1* ^ K ^ K I 

where v(n) is the perturbed eastward wind at trajectory point n, v(n) is 
the mean east wind at point n, v R (n-l) is the east wind perturbation 
from the previous point, P R (n) is either the large or small scale 
density perturbation calculated from (2). The coefficients I,J, and K 
are given by: 


I = 

K 


[o- (n)/<r (n-1)] [ (R - R R R )] 

V V VK £K V2£2 Vipi 

(1-R 2 R 2 ) 
pn vipi 


( 21 ) 




( 22 ) 


K t « [*>) - I K % 2 (n-1 ) - JV(.) 


' 2 Vi:VvlplV n ' 1 , V (n)1 


1/2 


(23) 


where R vr = R uk , and R y2 p 2 and R y are taken from tables. 

In equation (5), for R^, there is one scale distance for vertical 
displacements, L , and another for horizontal displacements, L . The 

* H 

small scale correlation, R^, is defined using another pair of scale 
distances. Similarly, there are vertical and horizontal scale distances 
for each of the temperature correlations. The large scale, horizontal 
scale distance, L , is a linear function of altitude: 

H 

L 

L = 900 + 6z (24) 

H 

L 

where z is the altitude and both z and L are in kilometers. The small 

H 

L 

scale, horizontal scale distance is given by 

L = 20 + 0.0125Z 2 (25) 

H 
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The vertical scale distances are functions of latitude as well as 


of altitude, specifically 


L y = [a + b(90 - [0. 22 + 0.00258z l,S ] 


(26) 


where <f> is the latitude in degrees, and a and b are empirical 

coefficients. The coefficients a and b are different for large and for 

small scale and also different for p and T. Thus there are different 

L , L , L , L . The horizontal scale distances are the same for 
v v v v 
pt ps TL TS 

density and temperature; L = L and L = L . The 

H H H H 

PS TS PL TL 

coefficients a and b are given in Table 1. 

Note that in (2) and (7) (and similarly for the small scale 
temperature and density perturbations) there is a gaussian random term 
which is added directly to the perturbation. Since gaussian random 
variables can take on arbitrarily large values, it is clearly possible 
that the perturbations can become unbounded. This could force the 
atmosphere attributes to become negative, which is physically 
unrealizable. In the GRAM code this modelling error is compensated by 
placing the following lower bounds on density, pressure and temperature 
perturbations: 


p’ = -0.9 p 

P’ = -0.9 P (27) 

T’ = -0. 9 T 

where, 

P’ =P L + P 5 

P’ = P + P (28) 

L S 

T* = T + T 

L S 

Thus, the density, temperature, and pressure are never allowed to fall 
below 1/10 of their respective mean values. Unfortunately, this does 
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not guarantee that pressure, density, and temperature will satisfy the 
perfect gas law. This problem is resolved in the current work by 
recalculating the perturbations whenever one is realized below 1/10 of 
its mean. Only the density and temperature are checked since the 
pressure is obtained by the perfect gas law from temperature and 
density. Either of these approaches will slightly skew the statistics 
by raising the mean and reducing the standard deviation. 

In the GRAM program as implemented, there can be discontinuities in 
the winds between 90 and 115 km altitude. In this range, the GRAM model 
fairs between the Groves section of the model, which is taken from 
tables, and the Jacchia section, which is analytic. Shown in figure 1 
are the magnitudes of north wind at three altitudes, over a range of 
latitudes and longitudes. These are mean wind magnitudes taken from a 
series of GRAM runs at constant altitude and latitude and varying 
longitude from 0 to 360°. At 95 km the wind velocity exhibits 
step changes with small changes in longitude. Just below this faired 
region, at 89 km, and above this region, at 116 km altitude, the winds 
vary continuously over the entire area. The winds at 95 km are 
representative of the winds at other altitudes in this region. Figure 2 
shows that the same phenomenon occurs in the eastward wind velocities. 
Because of this discontinuity and the fact that the latest modification 
to GRAM deletes winds above 125 km, winds are not used above 90km in the 
present implementation. 


IMPLEMENTATION OF SIMULATION IN POST 

This section describes the software implementation of the preceding 
model. The three main points to be considered are: how the theory is 
implemented into the existing POST program, the code changes necessary, 
and the results of these changes. 

The atmosphere model does not employ the standard POST atmosphere 
setup. A Runge-Kutta integrator must evaluate the state at several 
times in each integration step, effectively taking a finite difference 
derivative across each integration step. Since the POST atmospheres are 
easily differentiable, this integration scheme is entirely appropriate. 
An atmosphere model which includes Gaussian random terms is not 
differentiable except in a least squares sense. Because of this, use of 
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an atmosphere model which varies randomly within an integration step 
will lead to unpredictable results and possibly divergence of the 
vehicle simulation. 

Note that for a stochastic system, trajectory optimization is not 
defined in the usual sense. None of the optimization options in POST 
will be meaningful while using the random atmosphere. 

POST is a very large and complicated program. Therefore, it is 
desirable to make as few changes as possible to the core code. Adding a 
random atmosphere function in the existing framework would involve 
rerouting of the flow of control in several places. 

In this implementation, the random atmosphere is called through 
POST’s closed-loop guidance routine, CLGM. This insures that the 
atmosphere is updated at every time step. Since the CLGM routine is 
only called once per time step, the atmospheric terms remain constant 
across each integration interval, assuring well-behaved state 
integration. Thus implemented, the random atmosphere requires only one 
change to POST’s flow of control. 

The required changes are of four types: addition of new variables, 
addition of new tables, accessing new subroutines, and required settings 
to run the random atmosphere. The POST manual, [5, pp. 2-19 through 
2-22] describes the addition of new tables and variables. The specific 
additions are discussed in appendix A. 

All of the new variables added were placed in a new common block, 
called RAND. The variables in the common block RAND are X0LD(4) and 
AR0LD(12), where XOLD contains the state information at the last updated 
trajectory point, latitude, longitude, altitude, and dot product of 
downrange and crossrange. AROLD contains the mean and standard 
deviations of air density, pressure, and temperature for both large and 
small scale and the realizations of random pressure, density, and 
temperature at the last trajectory point. AROLD also contains the most 
recent realized perturbations in density (p ,p ) and temperature 

<\-V- 

Six new tables were added to the program. One each for large scale 
standard deviations and small scale standard deviations of pressure, 
density, and temperature. These tables may be functions of one, two, or 
three variables. The standard deviations used are from the GRAM-86 
model [2]. If other tables of standard deviations were available, they 
could be substituted in the POST NAMELIST. The GRAM data give each as 
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a function of altitude and latitude. 

A mean atmosphere is supplied as a set of three tables [4, pp. 

6. c-1 through 6.d-9]. This table need not have a regular grid; that is, 
it can have a larger number of data points in the most critical areas in 
the atmosphere and have relatively few data points in the less critical 
areas. The atmosphere table can be one, two, or three dimensional. 

There are a few settings for the POST NAMELIST that need to be made 
for the random atmosphere case. 

1. NPC(5)=0: This is the "no atmosphere" option. An atmosphere table 

is still input by the user and read as if NPC(5)=1 (the user supplied 
atmosphere) but the flag should be set to 0 to prevent the main POST 
atmosphere routine from overwriting the random values. [4,pg. 6. a. 4-1] 

2. NPC(12)*0: This insures that the downrange and crossrange will be 

calculated, therefore the distance traveled in each time step is known 
and the atmosphere correlation is calculated correctly. Typically, 
NPC(12)=2. [4, pg. 6. a. 19-1] 

3. MDL=8: This directs the integration to continue until the event 
criterion is met or exceeded. Thus the program will not iterate to find 
the exact time a criterion was met, but will stop at the next integer 
time step. [4,pp. 6. a. 6-1 through 6-4] 

4. IGUID(14)=2: This invokes the closed loop guidance routine which 
contains the call to the random atmosphere routine. [4,pp. 6.b.7-l and 
7-2] 

5. NPC(2)=1: Uses the fourth order Runge-Kutta integration. (4,pg. 

6. a. 15-3] 

Three subroutines were added to the code. RNDATM is accessed 
through the existing routine CLGM as discussed previously. The others, 
RNDATI and XKR, are called by RNDATM. RNDATM is the random atmosphere 
generator. RNDATI is an initialization routine called only once from 
RNDATM. RNDATI gives values for the "previous" random point needed to 
start the recursive algorithm. The Gaussian random number generating 
subroutine, XKR, takes approximately two realizations from a uniformly 
distributed sequence and transforms them into a single realization from 
normal distribution [6]. 

On the first pass, RNDATM calls RNDATI and returns. RNDATI 
establishes initial values needed for the iterative formula used in 
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RNDATM. On successive calls, RNDATM computes pressure, density, 
temperature and speed of sound as described above. 

The result of the RNDATM subroutine is an atmosphere profile which 
varies from the mean atmosphere in a random fashion according to the 
GRAM statistics. 


RESULTS 

The traditional method for doing Monte Carlo simulations through a 
random atmosphere involves precomputing a series of random atmosphere 
altitude profiles and using these profiles in the vehicle simulation. 
This procedure is illustrated in Figure 3. First the vehicle 
simulation is run using a deterministic atmosphere model. This produces 
a ’nominal’ trajectory which is edited to be monotonic in altitude and 
to insure the data points are properly spaced. This nominal trajectory 
is input into GRAM repeatedly to produce a series of atmosphere profiles 
which are dependent on altitude. These profiles are then reformatted to 
be used in the vehicle simulation. Reference 10 provides an excellent 
application of this approach. 

The atmosphere model described in previous sections of this report 

propagates an atmosphere state which depends continuously on the vehicle 

state and a vehicle state which is a function of the atmosphere state. 

Use of this model in Monte Carlo simulations will be referred to as the 

copropagated approach. Figure 4 shows the general flow of the 

copropagated approach. The POST plant model determines the vehicle 

state as a function of time while the atmosphere model determines the 

atmosphere state as a function of location, r. In the figure, l/z t 

represents the backward time-shift operator, and 1/z represents the 

r 

backward spatial-shift operator. The vehicle plant and the atmosphere 
model are connected by determining the spatial trajectory step, Ar, as a 
function of the time step. The atmosphere model is also driven by a 
white noise term, w(r). 

To demonstrate the copropagated atmosphere trajectory 
implementation, an aeroassisted orbital transfer maneuver is simulated. 
The purpose of such a maneuver is to utilize aerodynamic forces and the 
energy losses associated with them to change a vehicle’s orbit, either 
in radii, inclination, or both. The guidance law incorporated here was 
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developed for the Aeroassisted Flight Experiment (AFE) and treats the 
trajectory as one made up of two phases [7], The first, or entry phase 
is considered to be an equilibrium glide in which the rate of altitude 
change, dh/dt, is held constant. The second, exit phase, begins near 
the lowest altitude of the trajectory and is controlled by a 
predictor-corrector algorithm which continually updates control variable 
values to ensure a desired exit velocity and flight path angle. 

Statistical information was obtained by performing Monte Carlo 
analysis of randomly varying atmospheres encountered during 100 
simulated AFE flights. The results that follow include both trajectory 
and point statistics. Point statistics at atmosphere exit are compiled 
to illustrate the accumulated effects of the random atmosphere 
variations propagating throughout each flight (Table 2). Trajectory 
statistics are the result of averaging atmospheric variables (pressure, 
density, temperature) at each time step over the 100 flights to yield an 
"average" trajectory which is truncated to the minimum sample trajectory 
time duration. For presentation, these data are normalized about 
variable values corresponding to identical AFE trajectories simulated in 
a ’76 Standard Atmosphere without using the copropagated 
Atmosphere/Trajectory capability. 

Figure 5 shows altitude histories for both the average copropagated 
Atmosphere/Trajectory case and the same AFE profile employing a ’76 
Standard atmosphere. The shift in altitude developing in the latter 
portion of the flight can be attributed to a bias in air densities found 
between the ’76 and GRAM models effectively allowing the copropagated 
runs to generate greater lift in a more dense atmosphere. Plotted as a 
function of velocity, we see in figure 6 the same average trajectory and 
’76 Standard trajectory. This is a good illustration of the different 
spacing of data points during both the entry and exit portions of the 
flights. We also notice here the near constant velocity entrance and 
exit bracketing a region of large velocity gradients at lower altitudes. 

Normalized values of pressure, density and temperature using the 
copropagated atmosphere are plotted as functions of altitude in figures 
7, 8, and 9. These atmospheric variables display smoother structure of 
mean value fluctuations during entry (122 km altitude) than on exit. 

This can be attributed to the wider spacing of data points in this 
region caused by the high entrance velocity. Any shift in means at a 
given altitude between the entrance and exit phases is partially a 
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function of the latitude dependence of the GRAM atmosphere model, as the 
sample flights involved traversing a range of both latitudes and 
longitudes. 

Variances of normalized atmospheric variables are plotted in 
figures 10, 11, and 12. On all plots variances start at zero, 
consistent with the output of GRAM which assigns only mean values to 
initial trajectory locations. When plotted against altitude, variances 
display a pattern indicating that towards the end of the AFE flight 
variance magnitude and fluctuation increase. Trends in variance 
altitude histories display the effects of spatial separation. It is 
clear that when large distances separate consecutive trajectory points 
the GRAM perturbation model yields a smooth transition from one point to 
the next. 

A comparison of point statistics pertaining to Monte Carlo analysis 
is presented in. Table 2. The copropagated atmosphere and trajectory 
determination method is compared to the precomputed atmosphere profile 
approach. The mean trajectories for the two approaches were nearly 
identical. The final statistics of the two approaches (energy, 
latitude, and longitude) are in agreement, but there is a discrepancy in 
maximum dynamic pressure and heating rate and their corresponding 
variances. This is because at the near-constant altitude portion of the 
trajectory, where maximum dynamic pressure and heating rate are expected 
to occur, much greater density variance is modeled when using the 
copropagated capability. When determining trajectory and atmospheric 
variables simultaneously, variances are calculated based on the 
magnitude of spatial separation. The precomputed approach bases 
variances solely on vertical spatial separation which, in the portion of 
the trajectory being discussed here, is very small. This difference in 
the two techniques is made apparent in figures 13 and 14. Density 
variance is reduced nearly to zero at the midpoint of the trajectory 
when using the precomputed method, while the trend in copropagated 
variances appears to be much more realistic. The effects of this 
difference can be seen carrying over into the resulting velocity 
variances derived from the two approaches. These figures represent data 
generated in a previous though similar experiment using 80 simulated 
trajectories in the Monte Carlo analysis [9]. 
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CONCLUSIONS 


An add-on software module for POST which simulates spacecraft 
trajectories through a stochastic atmosphere based on the GRAM 
atmosphere model has been developed and demonstrated. This package 
gives realistic random fluctuations in density, pressure, and 
temperature while preserving the perfect gas relation. 

The software was written in such a way as to require only minor 
changes to the host program, POST. A complete listing of the changes to 
POST as well as the listings of the subroutines added are included in 
the appendices. 

This software package was exercised in simulating an aeroassisted 
orbit transfer problem and the results were compared to a simulation in 
which precomputed random atmosphere profiles were used. A comparison of 
results demonstrated significant differences in the statistics of key 
trajectory parameters between the two simulation approaches. This 
strongly suggests that the copropagated approach is to be preferred for 
Monte Carlo simulation, since it is capable of representing atmosphere 
variation in horizontal flight which the precomputed approach failed to 
exhibit. 
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APPENDIX A: Addition of Common Blocks and Tables to POST 


This appendix describes the changes made to the POST code to allow 
the use of the random atmosphere simulation. The procedure followed is 
as described in [5] except as required for the VAX operating system. 

The common block RAND was added as follows: 

1. The common block was referenced in subroutine DATA in file EXEC. FOR. 

In DATA a data statement was added that gave a default value for each 
variable in the common block. 

DATA X0LD/4*0. / 

DATA AR0LD/16*0. / 

2. In subroutine DICT, in file EXEC. FOR, each variable was assigned a 
Hollerith value. In subroutine DICT the name of the common block is 
different from that in DATA. In DICT the common was called RANDD but 
the variables were the same and in the same order. When the common was 
referenced in the execution of the code, it was referenced as RAND. 

POST has a general table facility which allows the addition of new 
tables without adding dimensioned arrays, hard-wired table arguments, 
table types, table dimensions, etc. All of the user input tables are 
placed in an array in blank common. Each table is referenced by a 
pointer. This pointer directs a general table lookup routine (GENTAB) 
to the proper location in the blank common. Each table also has a 
numeric and a mnemonic multiplier associated with it. Thus, to add a new 
table one needs only add a pointer and two multipliers to the comdeck 
MOTBL. 

1. A new common block, MOTBLO, was added to comdeck MOTBL. It contains 
the pointers for six tables: SDPRST, for small scale pressure standard 
deviation, SDPRLT, for large scale pressure standard deviation, SDRHST, 
SDRHLT, SDTMST, and SDTMLT, for density and temperature similarly. It 
also contains the multipliers for each of these tables: SDPRSM(2), etc. 
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2. An exact duplicate of this common, called MOTBLDO, was added to 
comdeck MOTBLD in file MOTBLD. INC. 

3. Data statements were added to subroutine DICT (which referenced 
common MOTBLDO) to give Hollerith values to the pointers and multipliers 
as: 

DATA SDPRST , SDPRSM/ 6HSDPRST , 6HSDPRSM , 5H=>0NE/ 

etc. 

4. In subroutine RTBLML (file RTBLML. NML) the table multipliers were 
added to the namelist. The data statement for NAMTBL was extended to 
include the six multiplier names: 

SDPRSM, SDPRLM, SDRHSM, SDRHLM, SDTMSM, SDTMLM. 

The data statement for ATRTBL was extended to make room for the new 
tables. 
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APPENDIX B: Source Listings 


The following is a source listing of the three subroutines added to 
POST for the random atmosphere. The first, RNDATM is called from the 
POST routine CLGM. RNDATM in turn calls the other two routines. The 
INCLUDE statements all access files used in the ordinary POST 
compilation, except for RAND. INC which is described in Appendix A. The 
scale height functions were read from the GRAM code from subroutine 
ERTRB. 
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SUBROUTINE RNDATM 
IMPLICIT REAL* 8 ( A-H , 0-Z ) 

INTEGER* 4 IX 

REAL *8 LVLRHO , LVSRHO , LVLTMP , LVSTMP , LVLWND , LVSWND 
REAL* 8 LHL , LHS 

** RNDATM *** 


NOTE: THE FORMULAS USED IN THIS SUBROUTINE ARE FOUND 
IN THE GRAM-86 REPORT BY JUSTUS, ET AL . 

MAKE SURE NPC(12)<>0 


INCLUDE ' POST3D_INCLUDE : MOTBL . INC/LIST ' 
INCLUDE' POST3D_INCLUDE: NPC . INC/LIST' 
INCLUDE ' POST3D_INCLUDE : GCRADC . INC/LIST ' 
INCLUDE 'POST3D_INCLUDE: MOTIC . INC/LIST' 
INCLUDE ' POST3D_INCLUDE : MOTVC . INC/LIST ' 
INCLUDE ' POST3D_INCLUDE : SERVC . INC/LIST ' 
INCLUDE 'POST3D_INCLUDE: RAND. INC/LIST' 
INCLUDE 'POST3D_INCLUDE: GUIDJC . INC/LI ST ' 
INCLUDE 'POST3D_INCLUDE: SPECAL . INC/LIST' 
INCLUDE ' POST3D_INCLUDE : AUXVC . INC/LIST' 
INCLUDE ' POST3D_INCLUDE: ALTITO. INC/LIST' 
INCLUDE 'POST3D_INCLUDE: TIME . INC/LIST ' 
INCLUDE ' POST3D_INCLUDE : VXIVC . INC/LIST' 
INCLUDE' POST3D_INCLUDE: VELOVC . INC/LIST' 
INCLUDE ' POST3D_INCLUDE : FUDGE . INC/LI ST ' 

***EMQ*** 

TRIAL COMMON 

COMMON/TRIAL/WORK ( 50 ) 


C 


C 


DIMENSION A( 3 ) ,B( 3 ) ,C( 3 ) 

EQUIVALENCE ( XOLD ( 1 ) , OLAT ) 
EQUIVALENCE ( XOLD ( 2 ) , OLONG ) 
EQUIVALENCE ( XOLD ( 4 ) , OLDZ ) 

EQUIVALENCE ( AROLD ( 1 ) , OLDPRS ) 
EQUIVALENCE ( AROLD ( 2 ) , OLDRHO ) 
EQUIVALENCE ( AROLD ( 3 ) , OLDTMP ) 
EQUIVALENCE ( AROLD ( 4 ) ,OSPRSS) 
EQUIVALENCE ( AROLD ( 5 ) , OSRHOS ) 
EQUIVALENCE ( AROLD ( 6 ) , OSTMPS ) 
EQUIVALENCE ( AROLD ( 7 ) , OGMPRS ) 
EQUIVALENCE ( AROLD ( 8 ) , OGMRHO ) 
EQUIVALENCE ( AROLD ( 9 ) , OGMTMP ) 
EQUIVALENCE ( AROLD (10), OSPRSL ) 
EQUIVALENCE ( AROLD ( 11 ) , OSRHOL ) 
EQUIVALENCE ( AROLD (12), OSTMPL ) 
EQUIVALENCE ( AROLD (13), ORHOPL ) 
EQUIVALENCE ( AROLD (14), OTMPPL ) 
EQUIVALENCE ( AROLD (15), ORHOPS ) 
EQUIVALENCE (AROLD (16) ,OTMPPS) 
EQUIVALENCE ( AROLD (17), OSUL ) 
EQUIVALENCE ( AROLD ( 18 ) ,OSUS) 
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EQUIVALENCE ( AROLD (19), OS VL ) 

EQUIVALENCE (AROLD( 20 ) , OSVS ) 

EQUIVALENCE (AROLD (21) ,OCURHL) 

EQUIVALENCE ( AROLD ( 22 ) , OCURHS ) 

EQUIVALENCE ( AROLD (23), OCVRHL ) 

EQUIVALENCE ( AROLD (24), OCVRHS ) 

EQUIVALENCE ( AROLD ( 25 ) ,OUPL) 

EQUIVALENCE ( AROLD (26), OUPS ) 

EQUIVALENCE ( AROLD ( 27 ) ,OVPL) 

EQUIVALENCE (AROLD (28) ,OVPS) 

C 

IF(NPC(5) .NE.O.AND.IGF(l) .EQ.O) RETURN 
if (npc(5) .ne.0)goto 200 
C 

IF( IGF( 1 ) . EQ. 0 )THEN 
CALL RNDATI 

ISEED»-INT(100.0*SECNDS(0.0) ) 

ISEED-INT( 1 . 0D8*RAN0 ( I SEED ) ) 

iseed«NINT( SPECI ( 8 ) ) 

rbar-0.0 

ictr-0 

RETURN 

END IF 

C JATM DETERMINES WHICH ATMOSPHERE TO USE AS THE MEAN 
JATM-1 

IF( JATM.GE. 2) GOTO 20 
C MEAN ATMOSPHERE INPUT IN TABLES 
GMPRS-EXP ( GENTAB ( PREST ) ) 

GMTMP-GENTAB ( ATEMT ) 

GMRHO-EXP ( GENTAB ( DENST ) ) 

PRES-GMPRS 

ATEM-GMTMP 

DENS-GMRHO 

C MEAN WINDS FROM TABLES 
GMU-GENTAB ( VWUT ) 

GMV-GENTAB ( VWVT ) 

VW(1)«GMU 
VW ( 2 ) **GMV 
VW ( 3 ) -GENTAB ( VWWT ) 

C VW(1)=0.0 

C VW( 2 )-0 . 0 

C VW( 3 )-0 . 0 

CS = DSQRT( ATMOSK( 1 ) *GMPRS ) 

GOTO 100 

20 IF( JATM.GE. 3 ) GOTO 30 

C 1962 STANDARD ATMOSPHERE AS MEAN 
CALL ATMOS 2 
GOTO 100 

30 IF( JATM.GE. 4) GOTO 40 

C 1963 PATRICK AFB ATMOSPHERE AS MEAN 
CALL ATMOS 3 
GOTO 100 

40 IF( JATM.GE. 5) GOTO 50 
C 1971 VANDENBURG AS MEAN 
CALL ATMOS 4 
GOTO 100 
50 CONTINUE 

C 1976 STANDARD ATMOSPHERE 
CALL ATMOS 5 
100 CONTINUE 
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IF ( SPEC I ( 9 ) . LT .0.0) GOTO 200 

C******************************************************************* 
C * * 

C * * 

I F ( GMPRS . LT . 0 . 0 ) GMPRS-0 . 0 

'C 


C GET NEW STD. DEV. FROM TABLE 
C 

SDTMPL-DSQRT ( GENTAB ( SDTMLT ) ) *GENTAB ( SDTMST ) 
SDPRSL-DSQRT ( GENTAB ( SDPRLT ) ) *GENTAB ( SDPRST ) 
SDRHOL-DSQRT ( GENTAB ( SDRHLT ) ) *GENTAB ( SDRHST ) 

C 

SDTMPS-DSQRT ( FPl-GENTAB ( SDTMLT ) ) *GENTAB ( SDTMST ) 
SDPRSS-DSQRT ( FPl-GENTAB ( SDPRLT ) ) *GENTAB ( SDPRST ) 
SDRHOS-DSQRT ( FPl-GENTAB ( SDRHLT ) ) * GENTAB ( SDRHST ) 

C 

C TYPE I CORRELATION 

C SDUL-DSQRT ( GENTAB ( SDULT ) ) *GENTAB ( SDUST ) 

C SDVL-DSQRT( GENTAB (SDVLT) ) *GENTAB ( SDVST ) 

C SDUS=DSQRT( FPl-GENTAB ( SDULT) ) *GENTAB( SDUST) 

C SDVS-DSQRT( FP1-GENTAB( SDVLT) )* GENTAB ( SDVST ) 

C 

C TYPE II CORRELATION 

SDUL-DSQRT ( GENTAB ( SDULT )*DABS( GENTAB (SDUST) ) ) 
SDVL-DSQRT ( GENTAB ( SDVLT )*DABS( GENTAB (SDVST) ) ) 
SDUS=DSQRT ( FPl-GENTAB ( SDULT ) *DABS ( GENTAB ( SDUST ) ) ) 
SDVS-DSQRT( FPl-GENTAB ( SDVLT ) *DABS ( GENTAB ( SDVST ) ) ) 


C 

C GET DENSITY-VELOCITY CORRELATIONS FROM TABLE 
C 


CURHOL-GENTAB ( CURHLT ) 
CURHOS-GENTAB ( CURHST ) 
CVRHOL-GENTAB ( CVRHLT ) 
CVRHOS-GENTAB ( CVRHST ) 




C * * 

C * * 


C 


I F ( GMRHO . LT . 0 . 0 ) GMRHO-0 . 0 


C 

C* 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 


** CALCULATE NEW RANDOM DENSITY 

CONVERT FEET TO KILOMETERS 
FTKM-12. 0/39370.0 

CONVERT NAUTICAL MILES TO KILOMETERS 
XNMIKM-6076 . 115*FTKM 

CALCULATE HORIZONTAL DISTANCE TRAVELED USING POSITION VECTOR 
DO 3 KJJ-1,3 

A( KJJ ) “XOLD( KJJ ) 

3 B( KJJ) -XI (KJJ) 

CALL VCROSS ( A, B , C ) 

CALL VUNI T ( A , A , XMAG ) 

CALL VUNIT(B,B,XMAG) 

CALL VUNI T(C,C, XMAG) 

IN ORDER TO AGREE WITH GRAM WE CALCULATE DISTANCE TRAVELED 
USING THE FLAT-EARTH THEORY 
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C PER « POLAR EARTH RADIUS IN KM 
C EER - EQUATORIAL EARTH RADIUS IN KM 
C EPS = ECCENTRICITY 
PER=6356 .7747 
EER-6378.160 

EPS«=( FPl-PER*PER/(EER*EER) ) 

RADEG-3 . 1415927D0/180 . 0D0 
DLON- ( LONG-OLONG ) *RADEG 
RGDLAT=GDLAT*RADEG 
DLAT- ( GDLAT-OLAT ) *RADEG 

DELX«GCRAD*DSQRT( DLAT**2+( COS ( RGDLAT ) *DLON ) * *2 ) 

I F ( IOFLAG.NE. 3 )DELX-DELX*FTKM 
DELZ=DABS ( OLDZ-ALTITO ) 

I F ( IOFLAG.NE. 3 )DELZ-DELZ*FTKM 

Z*»ALTITO*FTKM 

PHEE-DABS(GCLAT) 

DPHI-( 90 . 0-PHEE ) **2 
FZ-0. 22+0. 00258*Z* *1.5 
IF(FZ.GT. 5.0D0)FZ-5.0D0 

LARGE-SCALE SCALE HEIGHT - DENSITY 
LVLRHO™ (20.7-1. 346E-3*DPHI ) *FZ 
LARGE-SCALE SCALE HEIGHT - TEMPERATURE 
LVLTMP-7 . 3*FZ 

LARGE-SCALE SCALE HEIGHT - WIND 

LVLWND- ( 31 . 2-3 . 503E-3*DPHI ) *FZ 
SMALL-SCALE SCALE HEIGHT - DENSITY 
LVSRHO- ( 11 . 0-2 . 102E-4*DPHI ) *FZ 
SMALL-SCALE SCALE HEIGHT - TEMPERATURE 
LVSTMP= (3.0+5. 146E-4*DPHI ) *FZ 
SMALL-SCALE SCALE HEIGHT - WINDS 

LVSWND- (6.2-3 . 615E-4*DPHI ) *FZ 
LARGE-SCALE HORIZONTAL SCALE 
LHL=900 . 0+6 . 0*Z 
SMALL-SCALE HORIZONTAL SCALE 
LHS-20 . 0+0. 0125*Z**2 
IF(LHS.GT.400.0D0)LHS-400.0D0 

TERMS NEEDED FOR CALCULATING PERTERBATIONS 

RRHOL-EXP ( -1 . 0 *DSQRT ( ( DELX/LHL ) * * 2+ ( DELZ/LVLRHO ) * * 2 ) ) 
RTEEL«EXP(-1.0*DSQRT( ( DELX/LHL ) **2+ ( DELZ/LVLTMP )* *2 ) ) 
RUL-EXP( -1 . 0*DSQRT( ( DELX/LHL ) **2+ ( DELZ/LVLWND ) **2 ) ) 

RRHOS=EXP ( -1 . 0*DSQRT( ( DELX/LHS ) **2+ ( DELZ/LVSRHO) **2 ) ) 
RTEES“EXP ( -1 . 0*DSQRT( ( DELX/LHS ) **2+ ( DELZ/LVSTMP ) **2 ) ) 
RUS-EXP ( -1 . 0*DSQRT( ( DELX/LHS ) * *2+ ( DELZ/LVSWND ) **2 ) ) 

RVL-RUL 

RVS*=RUS 

AL-RRHOL*SDRHOL/OSRHOL 
BL-SDRHOL*DSQRT( FPl-RRHOL**2 ) 

AS=RRHOS*SDRHOS/OSRHOS 
BS-SDRHOS*DSQRT( FPl-RRHOS**2 ) 

10 CONTINUE 

*** CALCULATE NEW RANDOM TEMPERATURE 

RTRH02L=>( SDPRSL**2-SDRHOL**2-SDTMPL**2 )/( 2 . *SDRHOL*SDTMPL ) 
RTRHOlL-(OSPRSL**2-OSRHOL**2-OSTMPL**2)/(2.*OSRHOL*OSTMPL) 
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c 

RTRH02S=( SDPRSS**2-SDRHOS**2-SDTMPS**2 )/( 2 . *SDRHOS*SDTMPS ) 
RTRHOlS- ( OSPRSS**2-OSRHOS**2-OSTMPS**2 )/( 2 . *OSRHOS*OSTMPS ) 

C 

CEEL- ( SDTHPL/OSTMPL ) * ( ( RTEEL-RRH0L*RTRH02L*RTRH0lL )/ 

1 ( 1 .-RRHOL**2*RTRHOlL**2 ) ) 

CEES- ( SDTMPS/OSTMPS ) * ( ( RTEES-RRH0S*RTRH02S*RTRH0lS )/ 

1 ( 1 .-RRHOS**2*RTRHOlS**2 ) ) 

C 

DEEL-(RTEEL*SDTMPL-CEEL*OSTMPL)/(AL*RTRHOlL*OSRHOL) 

DEES® ( RTEES*SDTMPS-CEES*OSTMPS )/( AS* RTRHOlS* OSRHOS ) 

25 CONTINUE 

C I SEED® I NT ( 1 . 0D8 *RAN0 ( I SEED ) ) 

EEEL=DSQRT( SDTMPL**2-CEEL**2*OSTMPL**2-DEEL**2*SDRHOL**2-2 . 
1 *CEEL*DEEL*RRHOL*RTRHOlL*OSTMPL*SDRHOL ) 

EEES=DSQRT( SDTMPS**2-CEES**2*OSTMPS**2-DEES**2*SDRHOS**2-2 • 
1 *CEES*DEES*RRHOS*RTRHOlS*OSTMPS*SDRHOS ) 
c ****compute new perterbations in same order used by GRAM**** 
c***TOMO*** FORCED UNBIASED DETERMINATION of GAUSSIAN ERRORS (FUDGE) 
RBAR- ( ( RBAR*ICTR)+XKR( ISEED) )/( ICTR+1 ) 

ICTR-ICTR+1 
RLILl-XKR( ISEED)-RBAR 
PRINT* , ' FROM WITHIN RNDATM , ICTR - ',ICTR 
C RLILl**XKR( ISEED) 

RHOPS=GMRHO* ( AS* ( ( ORHOPS )/OGMRHO)+BS*RLILl ) 
c***TOMO*** 

RBAR=( ( RBAR* ICTR )+XKR( ISEED) )/( ICTR+1) 

ICTR- ICTR+1 
RLIL2=XKR( ISEED ) -RBAR 
C RLIL2=XKR( ISEED) 

TEMPPS-GMTMP* ( CEES* ( ( OTMPPS )/OGMTMP ) +DEES* 

1 ( (RHOPS)/GMRHO)+EEES*RLIL2) 
c 

**** calls for small scale U and V to keep number of random 
calls the same as GRAM***** 

***TOMO*** 

RBAR- ( (RBAR*ICTR)+XKR( ISEED) )/( ICTR+1 ) 

ICTR-ICTR+1 
RLIL3-XKR( ISEED)-RBAR 
RLIL3-XKR( ISEED) 

* * *TOMO* * * 

RBAR- ( (RBAR*ICTR)+XKR( ISEED) )/( ICTR+1) 

ICTR-ICTR+1 
RLIL4»XRR( I SEED) -RBAR 
RLIL4— XKR( ISEED) 

***TOMO*** 

RBAR-( (RBAR*ICTR)+XKR( ISEED) )/( ICTR+1) 

ICTR-ICTR+1 
RLIL5=XKR( ISEED ) -RBAR 
C RLIL5-XKR( ISEED) 

RHOPL-GMRHO* ( AL* ( ( ORHOPL )/OGMRHO ) +BL*RLIL5 ) 

C* * *TOMO* * * 

RBAR- ( (RBAR*ICTR)+XKR( ISEED) )/( ICTR+1) 

ICTR-ICTR+1 
RLIL6-XKR( ISEED)-RBAR 
C RLIL6=XKR( ISEED) 

TEMPPL-GMTMP* ( CEEL* ( ( OTMPPL ) /OGMTMP ) +DEEL* 

1 ( (RHOPL)/GMRHO)+EEEL*RLIL6 ) 
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c 

c **** calls for large scale U and V to keep number of random 
c calls the same as GRAM***** 

c 

c***tomo*** 

RBAR® ( ( RBAR* I CTR ) +XKR ( ISEED) )/( ICTR+1 ) 

ICTR=ICTR+1 
RLIL7»XKR( ISEED) -RBAR 
RLIL7-XKR( ISEED) 

***TOMO*** 

RBAR® ( ( RBAR*ICTR)+XKR( ISEED) )/( ICTR+1 ) 

ICTR-ICTR+1 
RLIL8-XKR( ISEED) -RBAR 
RLIL8-XKR( ISEED) 

WRITE ( 98,67)rlill,rlil2,rlil3,rlil4,rlil5,rlil6,rlil7,rlil8 
1 SDRHOS 

67 FORMAT (8F1 0.6) 

DENS-RHOPS+RHOPL+GMRHO 

DENS-GMRHO+RHOPS 

RHOLIM»0.1*GMRHO 

IF( DENS . LT. RHOLIM ) GOTO 25 

ATEM=TEMPPS+TEMPPL+GMTMP 

ATEM®GMTMP+TEMPPS 

TEMLIM=0.1*GMTMP 

I F ( ATEM . LT . TEMLIM ) GOTO 25 

*** CALCULATE NEW SPEED OF SOUND *** 

CS-DSQRT( ATMOSK( 1 ) *ATEM ) 

*** CALCULATE NEW RANDOM PRESSURE 

PRESPL®GMPRS* ( ( RHOPL/GMRHO ) + ( TEMPPL/GMTMP ) ) 

PRESPS=GMPRS * ( ( RHOPS/GMRHO ) + ( TEMPPS/GMTMP ) ) 

PRES® PRES PL+GMPRS+ PRES PS 
PRES-GMPRS+PRESPS 

*** CALCULATE RANDOM WIND PERTURBATIONS *** 

FL- ( SDUL/OSUL ) * ( RUL-RRHOL*CURHOL*OCURHL )/ 

1 ( FPl-RRHOL**2*OCURHL**2 ) 

FS® ( SDUS/OSUS ) * ( RUS-RRHOS*CURHOS*OCURHS )/ 

1 ( FPl-RRHOS**2*OCURHS**2 ) 

GL- ( RUL*SDUL-FL*OSUL )/( RRHOL*OCURHL*SDRHOL ) 

GS- ( RUS*SDUS-FS*OSUS )/( RRHOS*OCURHS*SDRHOS ) 

C 

HL-DSQRT( SDUL**2-FL*FL*OSUL**2-GL*GL*SDRHOL**2 
1 -FP2*FL*GL*RRHOL*OCURHL*SDRHOL*OSUL) 

HS-DSQRT( SDUS**2-FS*FS*OSUS**2-GS*GS*SDRHOS**2 
1 - FP2 * FS*GS*RRHOS*OCURHS* SDRHOS *OSUS ) 

C 

UPL-FL*OUPL+GL*RHOPL+HL*RLIL7 

UPS-FS*OUPS+GS*RHOPS+HS*RLIL3 

C 

XIL® ( SDVL/OSVL ) * ( RVL-RRHOL*CVRHOL*OCVRHL )/ 

1 ( FPl-RRHOL**2*OCVRHL**2 ) 

XIS® ( SDVS/OSVS ) * ( RVS-RRHOS*CVRHOS*OCVRHS )/ 

1 ( FPl-RRHOS**2*OCVRHS**2 ) 

C 
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XJL-(RVL*SDVL-XIL*OSVL)/(RRHOL*OCVRHL*SDRHOL) 

XJS- ( RVS*SDVS-XIS*OSVS )/( RRHOS*OCVRHS*SDRHOS ) 

XKL-DSQRT(SDVL**2-XIL*XIL*OSVL**2-XJL*XJL*SDRHOL**2 
1 -FP2*XI L*XJL*RRHOL*OCVRHL*SDRHOL*OSVL ) 

XKS«DSQRT( SDVS* * 2-XIS*XI S*OSVS* * 2-XJS*XJS*SDRH0S* * 2 
1 -FP2*XIS*XJS*RRH0S*0CVRHS*SDRH0S*0SVS ) 

VPL=XIL*0VPL+XJL*RH0PL+XKL*RLIL8 

VPS-XIS*0VPS+XJS*RH0PS+XKS*RLIL4 

WRITE ( 98 , * ) SDUS , OSUS , RUS , CURHOS , OCVRHS , SDVS , OSVS , CVRHOS 
1 , OCVRHS 

WRITE ( 98, *)FS,GS,HS,XIS,XJS,XKS 

WRITE ( 98 , * ) SDUL , OSUL , RUL , CURHOL , OCVRHL , SDVL , OSVL , CVRHOL 
1 , OCVRHL 

WRITE (98, * )FL,GL, HL,XIL,XJL,XKL 
VW( 1 ) -GMU* ( 1+UPL+UPS ) 

VW( 2 ) -GMV* ( 1+VPL+VPS ) 

** UPDATE 'OLD' VARIABLES 

XOLD ( 1 ) -GDLAT 
XOLD ( 2 ) *»LONG 
XOLD( 3 ) -XI ( 3 ) 

OLDZ-ALTITO 
AROLD ( 1 ) **PRES 
AROLD ( 2 ) -DENS 
AROLD ( 3 ) -ATEM 
AROLD( 4 )-SDPRSS 
AROLD ( 5 ) -SDRHOS 
AROLD ( 6 ) -SDTMPS 
AROLD ( 7 ) -GMPRS 
AROLD ( 8 ) -GMRHO 
AROLD ( 9 ) -GMTMP 
AROLD (10) — SDPRSL 
AROLD(ll)-SDRHOL 
AHOLD (12) -SDTMPL 
AROLD (13) -RHOPL 
AROLD (14) -TEMPPL 
AROLD (15) -RHOPS 
AROLD ( 16 )-TEMPPS 
AROLD ( 17 )-SDUL 
AROLD( 18 ) »SDUS 
AROLD (19) »*SDVL 
AROLD ( 20 ) -SDVS 
AROLD (21) -CURHOL 
AROLD( 22 )-CURHOS 
AROLD ( 23) -CVRHOL 
AROLD ( 24 ) -CVRHOS 
AROLD ( 25 ) -UPL 
AROLD( 26 ) -UPS 
AROLD ( 27 )-VPL 
AROLD ( 28 ) -VPS 
200 ITTT-NINT(TIME) 

SAVPRES— PRES 
SAVDENS-DENS 
SAVATEM-ATEM 
CALL ATMOS 5 
PRES76-PRES 
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ATEM76-ATEM 

DENS76=DENS 

PRES=SAVPRES 

DENS-SAVDENS 

ATEM=SAVATEM 

WLONG=360 . ODO-LONG 

WRI TE ( 9 9 , * ) I TTT , ALTI TO , GCLAT , WLONG , PRES , DENS , ATEM , VA , VW , 
1 PRES 7 6 , DENS 7 6 , ATEM76 

201 CONTINUE 
RETURN 
END 

SUBROUTINE RNDATI 
IMPLICIT REAL* 8 ( A-H , O-Z ) 

REAL* 8 LVLRHO , LVSRHO , LVLTMP , LVSTMP , LVLWND , LVSWND 
REAL* 8 LHL , LHS 
C 

C THIS INITIALIZES THE VARIABLES USED IN THE RANDOM 
C ATMOSPHERE MODEL 

INCLUDE ' POST3D_INCLUDE : RAND . INC/LIST ' 

INCLUDE ' POST3D_INCLUDE : MOTVC . INC/LIST ' 

INCLUDE' POST3D_INCLUDE: SPECAL . INC/LIST' 

INCLUDE ' POST3D_INCLUDE : MOTI C . INC/LIST ' 

INCLUDE ' POST3D_INCLUDE : AUXVC . INC/LIST' 

INCLUDE ' POST3D_INCLUDE : ALTITO. INC/LIST' 

INCLUDE ' POST3D_INCLUDE : MOTBL . INC/LIST ' 

INCLUDE ' POST3D_INCLUDE : VXIVC . INC/LIST ' 

INCLUDE ' POST3D_INCLUDE rGCRADC . INC/LIST' 

INCLUDE ' POST3D_INCLUDE : SERVC . INC/LIST ' 

INCLUDE ' POST3D_INCLUDE: FUDGE . INC/LIST' 

INCLUDE ' POST3D_INCLUDE : GUIDJC . INC/LIST ' 

C 

C XOLD ( 1 ) = OLD GEODETIC ALTITUDE 
C XOLD( 2 ) = OLD LONGITUDE 
C XOLD( 3 ) = INERTIAL POSITION VECTOR 
C XOLD ( 4 ) = ALTITUDE 
C 

C AROLD(l)- OLD PRESSURE 
C AROLD( 2 ) ■ OLD DENSITY 
C AROLD( 3 ) = OLD TEMPERATURE 

C AROLD( 4 ) = OLD PRESSURE STANDARD DEVIATION - SMALL SCALE 
C AROLD ( 5 ) = OLD DENSITY STANDARD DEVIATION - S.S. 

C AROLD ( 6 ) = OLD TEMPERATURE STANDARD DEVIATION - S.S. 

C AROLD ( 7 ) = OLD MEAN PRESSURE 
C AROLD ( 8 ) = OLD MEAN DENSITY 
C AROLD ( 9 ) “ OLD MEAN TEMPERATURE 

C AROLD ( 10)= OLD PRESSURE STANDARD DEVIATION - LARGE SCALE 
C AROLD (11)= OLD DENSITY STANDARD DEVIATION - L.S. 

C AROLD (12)= OLD TEMPERATURE STANDARD DEVIATION - L.S. 

C AROLD (13)= OLD L.S. DENSITY PERTURBATION 
C AROLD (14)= OLD L.S. TEMPERATURE PERTURBATION 
C AROLD (15)= OLD S.S. DENSITY PERTURBATION 
C AROLD (16)= OLD S.S. TEMPERATURE PERTURBATION 
C AROLD (17)= OLD L.S. NORTHERLY WIND STANDARD DEV. 

C AROLD ( 18 ) = OLD S.S. NORTHERLY WIND STANDARD DEV. 

C AROLD ( 19 ) = OLD L.S. EASTERLY WIND STANDARD DEV. 

C AROLD (20)= OLD S.S. EASTERLY WIND STANDARD DEV. 

C AROLD (21)= OLD L.S. DENSITY-NORTHERLY WIND CORRELATION 
C AROLD ( 22 ) = OLD S.S. DENSITY-NORTHERLY WIND CORRELATION 
C AROLD ( 23 )= OLD L.S. DENSITY-EASTERLY WIND CORRELATION 
C AROLD ( 24 ) = OLD S.S. DENSITY-EASTERLY WIND CORRELATION 
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C AROLD (25)- OLD L.S. NORTHERLY WIND PERTERBATION 
C AROLD( 26 ) - OLD S.S. NORTHERLY WIND PERTERBATION 
C AROLD ( 27 ) = OLD L.S. EASTERLY WIND PERTERBATION 
C AROLD (28)- OLD S.S. EASTERLY WIND PERTERBATION 
C 

C *** INITIALIZE MEAN VALUES *** 

C 

AROLD ( 7 ) -EXP ( GENTAB ( PREST ) ) 

AROLD ( 8 ) -EXP ( GENTAB ( DENST ) ) 

AROLD ( 9 ) -GENTAB ( ATEMT ) 

C 

IF(SPECI(9) .LT.0.0 )GOTO 200 
** INITIALIZE STANDARD ERRORS ** 


AROLD ( 4 ) -DSQRT ( FPl-GENTAB ( SDPRLT ) ) * GENTAB ( SDPRST ) 

AROLD ( 5 ) -DSQRT ( FPl-GENTAB ( SDRHLT ) ) *GENTAB ( SDRHST ) 

AROLD ( 6 ) -DSQRT ( FPl-GENTAB ( SDTMLT ) ) *GENTAB ( SDTMST ) 

AROLD (10) -DSQRT ( GENTAB ( SDPRLT ) ) * GENTAB ( SDPRST ) 

AROLD (11) -DSQRT ( GENTAB ( SDRHLT ) ) *GENTAB ( SDRHST ) 

AROLD (12) -DSQRT ( GENTAB ( SDTMLT ) ) *GENTAB ( SDTMST ) 

TYPE I CORRELATION 

AROLD (17) -DSQRT ( GENTAB ( SDULT ) ) *GENTAB ( SDUST ) 
i AROLD (18) -DSQRT ( FPl-GENTAB ( SDULT ) ) *GENTAB ( SDUST ) 

AROLD (19) -DSQRT ( GENTAB ( SDVLT ) ) * GENTAB ( SDVST ) 

AROLD ( 20 ) -DSQRT ( FPl-GENTAB ( SDVLT ) ) *GENTAB ( SDVST ) 

TYPE II CORRELATION 

AROLD (17) -DSQRT ( GENTAB ( SDULT ) *DABS ( GENTAB ( SDUST ) ) ) 
AROLD (18) -DSQRT ( FPl-GENTAB ( SDULT ) *DABS ( GENTAB ( SDUST ) ) ) 
AROLD ( 19 ) -DSQRT ( GENTAB ( SDVLT ) *DABS ( GENTAB ( SDVST ) ) ) 
AROLD (20) -DSQRT ( FPl-GENTAB ( SDVLT ) *DABS ( GENTAB ( SDVST ) ) ) 

AROLD ( 21 ) -GENTAB (CURHLT) 

AROLD (22) -GENTAB ( CURHST ) 

AROLD (23) -GENTAB ( CVRHLT ) 

AROLD ( 2 4 ) -GENTAB ( CVRHST ) 

*** LET INITIAL REALIZATIONS 

OF RANDOM COMPONENTS GIVE MEAN PLUS SQRT OF SUM OF SQUARES 
OF SMALL AND LARGE S.D. WITH GAUSSIAN RANDOM NUMBER.*** 


11 CONTINUE 

C PPRIME-DSQRT( AROLD( 10 ) **2+AROLD( 4 ) **2 ) *XKR( ISEED) 

PPRIME-0 . 0D0 

I F ( PPRIME. LT.-l . 0D0 )GOTO 11 
' AROLD( 1 )-AROLD( 7 ) * ( 1+PPRIME ) 

C 

22 CONTINUE 

C DPRIME-DSQRT( AROLD( 11 ) **2+AROLD( 5 ) **2 ) *XKR( ISEED) 

DPRIME-0 . 0D0 

I F ( DPRIME . LT . -1 . 0D0 ) GOTO 22 
AROLD ( 2 ) -AROLD ( 8 ) * ( 1+DPRIME ) 

C 

33 CONTINUE 

C TPRIME-DSQRT( AROLD ( 12 ) **2+AROLD( 6 ) **2 ) *XKR( ISEED) 
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TPRIME-0.0D0 

I F ( TPRIME . LT . -1 . 0D0 ) GOTO 33 
AROLD ( 3 ) - AROLD ( 9 ) * ( 1+TPRIME ) 

ASSIGN OLD PERTURBATIONS BASED ON RELATIVE SIZE OF S.D. 

SIGP-AROLD( 4 )+AROLD( 10 ) 

S I GD=AROLD ( 5 ) +AROLD ( 11 ) 

S I GT-AROLD ( 6 ) + AROLD (12) 

AROLD ( 1 3 ) « ( AROLD ( 1 1 ) /S I GD ) *DPRI ME * AROLD ( 8 ) 

AROLD ( 14 )»( AROLD ( 12 )/SIGT) *TPRIME*AROLD( 9 ) 

AROLD ( 1 5 ) « ( AROLD ( 5 ) /S I GD ) * DPRI ME * AROLD ( 8 ) 

AROLD ( 16 )«( AROLD ( 6 )/SIGT) * TPRIME* AROLD ( 9 ) 

AROLD ( 25 )"0 . 0D0 
AROLD ( 26 )®0 . 0D0 
AROLD ( 27 )=0.0D0 
AROLD ( 28 )“0 . 0D0 

SKIP TO HERE IF NOT USING PERTURBATTIONS 
200 CONTINUE 

XOLD ( 1 ) =GDLAT 
XOLD( 2 ) «LONG 
XOLD ( 3 ) =XI ( 3 ) 

XOLD( 4 )-ALTITO 

PRES=AROLD( 7 ) 

DENS-AROLD(8) 

ATEM«*AROLD( 9 ) 

CS-DSQRT ( ATMOSK ( 1 ) * AROLD ( 9 ) ) 

DON'T WNAT TO RE-INI TIALIZE 
IGF ( 1 ) “1 
RETURN 
END 

REAL*8 FUNCTION XKR(ISEED) 

GAUSSIAN RANDOM NUMBER GENERATOR 

THIS ROUTINE FOLLOWS THE ARTICLE BY KINDERMAN AND RAMAGE 
"COMPUTER GENERATION OF NORMAL RANDOM VARIABLES" 

IN THE JOURNAL OF THE AMERICAN STATISTICAL ASSOCIATION 
VOLUME 71, NUMBER 356, PP. 893-896. 

CODED BY ERIC QUEEN 7/30/88 

IMPLICIT REAL*8 ( A-H , O-Z ) 

REAL* 8 RANDOM 4 

PHI ( X)“DEXP( -0 . 5D0*X*X )/DSQRT( 2 . 0D0*3 .141 592654 DO) 
F(X)=PHI(X)-.180025191068563D0*( XI -DABS ( X) ) 

XI -2.216035867166471D0 
C BEGIN STEP 1 

U**DBLE ( RAN0 ( ISEED ) ) 

I F( U . LT .0.884070402298758 DO) THEN 
V=DBLE(RAN0( ISEED) ) 

XKR=XI* ( 1 . 131131635444180D0*U+V-1 . 0D0 ) 
RETURN 

ENDIF 

C BEGIN STEP 2 

IF ( U . LT . .973310954173898D0) THEN 
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c 


c 

c 

80 


C 


70 


C 


50 


C 

30 


BEGIN STEP 4 

I F ( U . LT . . 958720824790463D0 )THEN 

BEGIN STEP 6 

I F ( U . LT . . 911312780288703D0 )THEN 

BEGIN STEP 8 
V-DBLE ( RANO ( I SEED) ) 
W-DBLE( RANO ( ISEED) ) 

Z-V-W 

VWMIN-DMINl(V,W) 
VWMAX-DMAXl ( V , W ) 

T=. 4797274 042224 41D0-. 5955071 3801 5940D0*VWMIN 
I F ( VWMAX . LE . 0 . 805577924423817D0 )GOTO 90 
TEMPI-. 053377 549506886D0*DABS( Z) 

I F ( F(T) .GE. TEMPI )GOTO 90 
GOTO 80 


ELSE 


ENDIF 


BEGIN STEP 7 
V-DBLE ( RANO ( I SEED ) ) 

W-DBLE ( RANO ( ISEED ) ) 

Z-V-W 

VWMIN-DMINl ( V f W) 

VWMAX-DMAXl ( V , W ) 

T-.479727404222441D0+1. 105473661022070D0*VWMIN 
I F ( VWMAX . LE . . 872834976671790D0 )GOTO 90 
TEMP2-. 049264 496373128D0*DABS(Z) 

I F ( F ( T ) . GE . TEMP2 ) GOTO 90 
GOTO 70 


ELSE 


BEGIN STEP 5 
V-DBLE ( RANO ( ISEED ) ) 

W-DBLE ( RANO ( ISEED ) ) 

Z-V-W 

VWMIN-DMINl ( V,W) 

VWMAX-DMAXl (V,W) 

T=XI- . 630834801921960D0*VWMIN 
I F ( VWMAX . LE . . 755591531667601D0 )GOTO 90 
TEMP 3- . 03424050375011lD0*DABS ( Z ) 

I F ( F ( T ) . GE . TEMP3 ) GOTO 90 
GOTO 50 


ENDIF 


ELSE 


ENDIF 


BEGIN STEP 3 

V-DBLE(RAN0( ISEED) ) 

W-DBLE ( RANO ( ISEED) ) 

TEMP5-0. 5D0*XI*XI 
T-TEMP5-DLOG(W) 

TEMP4— V**2*T 
IF ( TEMP 4 . GT . TEMP 5 ) THEN 
GOTO 30 


ELSE 

I F ( U . LT . . 986655477086949D0 )THEN 
XKR-DSQRT( 2*T) 

RETURN 

ELSE 

XKR— DSQRT(2*T) 

RETURN 

ENDIF 

ENDIF 


29 



c 


BEGIN STEP 9 


90 CONTINUE 

IF ( Z . LT . 0 . 0D0 ) THEN 
XKR-T 
RETURN 

ELSE 

XKR=-T 

RETURN 

ENDIF 

RETURN 

END 

FUNCTION RANO(IDUM) 

DIMENSION V( 97 ) 

DATA IFF /0/ 

I F ( IDUM . LT . 0 . OR .IFF. EQ . 0 ) THEN 
IFF=1 

ISEED=ABS ( IDUM ) 

IDUM-1 

DO 11 J-1,97 
DUM=RAN ( I SEED) 

11 CONTINUE 

DO 12 J=1 , 97 

V( J )=RAN( I SEED) 

12 CONTINUE 
Y=RAN( ISEED) 

ENDIF 

J=1+INT ( 97 . * Y ) 

I F ( J . GT . 9 7 . OR . J . LT . 1 ) PAUSE 
Y=V( J) 

RANO-Y 

V(J)=RAN( ISEED) 

RETURN 

END 
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Table 1. Coefficients of Scale Distances 


Parameter 

a 

b 

L 

V 

pL 

20.7 

-1. 346x1 0" 3 

L 

V 

P s 

11.0 

-2. 102xl0 -4 

L 

V 

TS 

3.0 

5. 146xl0 -4 

L 

V 

TL 

7.3 

0.0 

-4 

L 

V 

us 

6.2 

3.615x10 

_ -3 

L 

V 

31.2 

-3.503x10 

(JL 



Table 2. Point 

Statistics from 

100 Sample AFE Trajectories 

Copropagated Precomputed 

/ m c & n f » / id c & n } . 

Ivarlance J (.variance; 


Maximum Dynamic Pressure 
(kg/m. sec ) 

1669.30 

10615.25 

1546.74 
6402. 82 

Maximum Heating Rate 
(btu/f t 2 sec) 

145.90 

19.69 

142.40 

11.33 

Heating Rate Integral 
(btu/f t 2 sec) 

20187.43 

101451.30 

19942. 69 
82547. 80 

2 2 

Final Energy (m /sec ) 

4232. 32 
73.97 

4231.67 

30.89 

Final Longitude (deg) 

33.0 

2.4 

32.41 

3.4 

Final Latitude (deg) 

24.32 
0. 16 

24.23 

0.24 

Minimum Altitude (km) 

74.25 

75.07 


0.035 0.049 
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Figure 3.- Traditional, precomputed Monte Carlo Simulation Approach 





Figure 4.- Copropagated Monte Carlo simulation approach 
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Fiqure 6 - Average Altitude Histories of AFE flight through 

both Copropagated and ‘76 Standard Atmospheres 




Figure 7.— Mean Pressures Normalized w.r.t. '76 Standard Atmosphere 
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Figure 8.— Mea 
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Figure 9. Mean Temperatures Normalized w.r.t. '76 Standard Atmosphere 



Figure 10. 
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Figure 11.— Variance of Normalized Densities 



Figure 12.— Variance of Normalized Temperatures 
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Figure 13.— Comparison of Density Variance Histories 

for Copropagated and Precomputed Approaches 
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Figure 14.— Comparison of Velocity Variance Histories 

for Copropagated and Precomputed Approaches 
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